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5_^ . Abstract 

^j^. Can we make Bayesian posterior MCMC sampling more efficient when faced with very large datasets? 

'^N ' We argue that computing the likelihood for iV datapoints twice in order to reach a single binary decision is 

O^ , computationally inefficient. We introduce an approximate Metropolis-Hastings rule based on a sequential 

^S) ' hypothesis test which allows us to accept or reject samples with high confidence using only a fraction of 
the data required for the exact MH rule. While this introduces an asymptotic bias, we show that this 

' ^' , bias can be controlled and is more than offset by a decrease in variance due to our ability to draw more 

V^ ' samples per unit of time. We show that the same idea can also be applied to Gibbs sampling in densely 

^J I connected graphs. 



1 Introduction 



s?J ' Markov chain Monte Carlo (MCMC) sampling has been the main workhorse of Bayesian computation ever 

Q"^ . since it's development in the 1950's. A canonical MCMC algorithm proposes samples from a distribution 

0^ ' q and then acc epts or rejects these pr oposals with a certain probability given by the Metropolis-Hastings 

OQ . (MH) formula Metropolis et al.l . Il953l | . The MH rule needs to examine the likelihood of all the data- items 



^i ' twice, once for the current state and another time the proposed state. When the number of data-cases is 

■^ . large this is an awful lot of computation for one bit of information, namely to accept or reject it. 

^^ ' In today's Big Data world, we need to rethink our Bayesian inference algorithms or risk becoming 

obsolete. Standard MCMC methods do not meet the Big Data challenge for the reason described above. 
Researchers have made some progress in terms of making MCMC more efficient, mostly by focussing on 
parallelization. Very few question the algorithm itself: is the standard MCMC paradigm really optimally 
efficient in achieving its goals? We claim its not. 
^ , Any method that includes computation as an essential ingredient should acknowledge that there is a finite 

^.' amount of time, T, to finish a calculation. An efficient MCMC algorithm should therefore decrease the "error" 

(properly defined) maximally in the given time T. For MCMC algorithms, there are two contributions to this 
error: bias and variance. Bias occurs because the chain needs to burn in during which it is sampling form 
the wrong distribution. Bias usually decreases fast, as evidenced by the fact that practitioners are willing 
to wait until the bias has (almost) completely vanished after which they discard these "burn-in samples". 
The second cause for error is sampling variance, which occurs because of the random nature of the sampling 
process. The retained samples after burn-in will reduce the variance roughly as 0{1/N). 

However, when given a finite time T it is not at all clear whether the strategy of retaining few unbiased 
samples and thus accepting an error dominated by variance is optimal. Perhaps, by decreasing the bias more 
slowly we could sample faster and thus reduce variance faster? In this paper we illustrate this effect by cutting 
the computational budget of the MH accept/reject step. We argue that many accept/reject decisions can 
be ta ken by examining only a fraction of the full dataset. To achieve that, we conduct sequential hypothesis 
tests [Govindaraiului l2004j to compute the probability of accepting or rejecting a given sample and find that 
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the majority of these decisions c an be made based on a small fraction of the data with high confidence. A 
similar method was developed in lSingh et alj 2012| . where the factors of a graphical model are sub-sampled 
to compute fixed width confidence intervals for the log-likelihood terms appearing in the MH test. 

Our "philosophy" runs deeper than the proposed algorithm here. We advocate MCMC algorithms with 
a "bias-knob" , allowing one to dial down the bias at a rate that optimally balances error due to bias and 
variance. We only know of one algorithm that would also adhere to this strategy: stochast ic gradient 
Langevin dynamics Welling and Tehl . 12011 1 and its successor stochastic gradient Fisher scoring Ahn et al.l . 
2012| . In their case the bias- knob was the stepsize. These algorithms do not have a MH step which resulted in 



occasional samples with extremely low probability. We show that our approximate MH step largely resolves 
this issue, still avoiding 0{N) computations per iteration. 

In the next section we introduce the MH algorithm and discuss its drawbacks. Then in section [31 we 
introduce the idea of approximate MCMC methods and the bias variance trade-off involved. We develop 
approximate MH tests in section |4] for Bayesian posterior sampling. We show that the same tests can be 
applied to Gibbs sampling [S] Finally, we show our experimental results in [5] and conclude in [T] 



2 The Metropolis-Hastings algorithm 

MCMC methods generate samples from a distribution 50(6*) by forward simulating a Markov chain designed 
to have stationary distribution So{9). A Marko v chain with a given st ationary distribution can be con- 
structed using the Metropolis- Hastings algorithm Metropolis et al.l . Il953j , which uses the following rule for 
transitioning from the current state 0f to the next state 9t+i: 

1. Draw a candidate state 9' from a proposal distribution q{9'\9t) 

2. Compute the accept probability: 

So{9')q{9t\9') 



P„ — niin 



1, 



S^{9t)q{9'\9t) 



(1) 



3. Draw u ~ Uniform[0, 1]. If u < Pa set 9t+i <— 9' , otherwise set 9t+i <— 9t. 

Following this transition rule ensures that the stationary distribution of the Markov chain is So{9). The 
samples from the Markov chain are usually used to estimate the expectation of a function f{9) with respect 
to So{9). To do this we collect T samples and approximate the expectation / — {f)so as I — ^ X]t=i f{^t)- 
Since the stationary distribution of the Markov chain is 5o, / is an unbiased estimator of /. 

The variance of / is V = &^{{f)so — ^f{(^t))'^], where the expectation is over multiple simulations of the 
Markov chain. It is well known that V ~ a'i ^ t/^, where a^f ^ is the variance of / with respect to Sq and 
T is the integrated auto-corr elation time, which is a measure of the interval between independent samples 
Gamerman and Lopesl . |2006| . Usually, it is quite difficult to design a chain that mixes fast and therefore, the 



auto-correlation time will be quite high. Also, for a lot of important problems, evaluating Sq{9) to compute 
the acceptance probability Pa in every step is so expensive that we can collect only a very small number 
of samples (T) in a given amount of computational time. Thus the variance of / can be prohibitively high, 
even though it is unbiased. 



3 Approximate MCMC and the Bias Variance Trade-ofF 

Ironically, the reason why MCMC methods are so slow is that they are designed to be unbiased. In fact, 
if we were to allow a sligh t bias in the stat i onary distribution, it is possible to design a Markov chain that 



can be simulated cheaply Welling and Tehl . l201l[|Ahn et al.l . l2012| . That is, to estimate / — (/) 



/5o' 



construct a Markov chain with stationary distribution S^ where e is a parameter that can be used to control 
the bias in the algorithm. Then / can be estimated as / = ^ St=i fi^t), computed using samples from 5^ 
instead of S^ . 

As e — >■ 0, Se approaches 5o (the distribution of interest) but it becomes expensive to simulate the Markov 
chain. Therefore, the bias in / is low, but the variance is high because we can collect only a small number 



of samples in a given amount of computational time. As e moves away from 0, it becomes cheap to simulate 
the Markov chain but the difference between S^ and the distribution of interest Sq grows. Therefore, / will 
have higher bias, but lower variance because we can collect a larger number of samples in the same amount 
of computational time. This is a classical bias-variance trade-off and can be studied using the risk of the 
estimator. 

The risk of / is defined as the mean squared error in /, i.e. R — K[{{f)s„ — ^/(6't))'^], where the 
expectation is taken over multiple simulations of the Markov chain. It is easy to show that the risk can be 
decomposed a.s R = B^ + V, where B = {f)s, - {f)so is the bias and V = E[((/)5, - 7/(6**))^] ~ ^},sj/'^ 
is the variance. 

The optimal setting of e that minimizes the risk depends on the amount of computational time available. 
If we have an infinite amount of computational time, we should set e to 0. Then there is no bias, and 
the variance can be brought down to because we can collect an infinite number of samples. This is the 
traditional MCMC setting. However, given a finite amount of computational time, this setting may not be 
optimal. It might be better to tolerate a small amount of bias in the stationary distribution if it allows us to 
reduce the variance quickly, either by making it cheaper to collect a large number of samples or by mixing 
faster. 

It is interesting to note that t wo recently proposed al gorithms follow this paradigm: Stochastic Gradi- 
ent Langev i n Dy namics (SGLD) [Welling and Tehl . I2OII and Stochastic Gradient Fisher Scoring (SGFS) 



Ahn et aU |2012| | . These algorithms are biased because they omit the required Metropolis-Hastings tests 



However, in both cases, a knob e (the step-size of the proposal distribution) is available to control the bias. 
As e — >■ 0, the acceptance probability Pa ^ 1 and the bias from not conducting MH tests disappears. How- 
ever, when e — > the chain mixes very slowly and the variance increases because the auto-correlation time 
T — >■ 00. As e is increased from 0, the auto-correlation, and therefore the variance, reduces. But, at the same 
time, the acceptance probability reduces and the bias from not conducting MH tests increases as well. 

In the next section, we will develop another class of approximate MCMC algorithms for the case where 
the target So is a Bayesian posterior distribution given a very large dataset. We acheive this by developing 
an approximate Metropolis-Hastings test, equipped with a knob for controlling the bias. Moreover, our 
algorithm has the advantage that it can be used with any proposal distribution. For example, our method 
allows approximate MCMC methods to be applied to problems where it is impossible to compute gradients 
(which is necessary to apply SGLD/SGFS). Or, we can even combine our method with SGLD/SGFS, to 
obtain the best of both worlds. 

4 Approximate Metropolis-Hastings Test for Bayesian Posterior 
Sampling 

An important method in the toolbox of Bayesian inference is posterior sampling. Given a dataset X^ con- 
sisting of N independent observations {xi, ...Xn} which we model using a distribution p{xi; 0) parameterized 
hy G R^, and a prior distribution p{0), the task is to generate samples from the posterior distribution 

Soie)^p{0)7rl,pix,;O). 

If the dataset has a billion datapoints, it becomes excruciatingly painful to compute So{.) in the MH test, 
which has to be conducted for each posterior sample that we generate. Spending 0{N) computation to get 
just 1 bit of information, i.e. whether to accept or reject a sample, is likely not the best use of computational 
resources. 

But, if we try to develop such accept/reject tests using only sub-samples of data, that satisfy detailed 
balance exactly with respect to the posterior distribution, w e will quickly see the no fre e lunch theorem 
kicking in. F or example, the p seudo marginal MCMC method Andrieu and Roberta . l2009l | and the method 



KicKmg m. ■f or example, tne p seuao margmai iviuiviu metnoa [Anarieu ana Koberta . l/UUalj ana tne metnoa 
developed by iLin et al.l 200Cll | provide a way to conduct exact accept/reject tests using unbiased estimators 



of the likelihood. However, unbiased estimators of the likelihood that c an be computed from mini-batches of 
data, such as the Poisson estimator or the Kennedy-Bhanot estimator iLin et al.l 2000l | have incredibly high 



variance for large datasets. Because of this, once we get a very high estimate of the likelihood, almost all 
proposed moves are rejected and the algorithm gets stuck. 

Therefore, we should be willing to tolerate some error in the stationary distribution if we want faster 
accept/reject tests. If we can offset this small bias by drawing a large number of samples cheaply and 



reducing the variance faster, we can establish a potentiahy large reduction in the risk. 

We will now show how to develop such approximate tests by reformulating the MH test as a statistical 
decision problem. It is easy to see that the original MH test is equivalent to the following procedure: Draw 
u ^ Uniform[0, 1] and accept the proposal 9' if the average difference fj, in the log-likelihoods of 9' and 9t is 
greater than a threshold /iq, i.e. compute 



Mo 



1 

N 



log 



Pi0t)q{9'\9t 



and 



pi9')q{9t\9') 
- ^ ;» where h = \ogp{xi;9') - \ogp{xi;9t) 



N 



(2) 



Then ii fi > fiQ, accept the proposal and set 9t+i ■'r- 9' . If /x < /xq, reject the proposal and set 9t+i <— 9t. This 
reformulation of the MH test makes it very easy to frame it as a statistical hypothesis test. Given /uq and 
a random sample {^i...Z„} drawn without replacement from the population {Zi.../jv}, can we decide whether 
the population mean // is greater than or less than the threshold ^q? The answer to this depends on the 
precision in the random sample. If the difference between the sample mean I and /xg is significantly greater 
than the standard deviation s of I, we can make the decision to accept or reject the proposal confidently. If 
not, we should draw more data to increase the precision of / (reduce s) till we have enough evidence to make 
a decision. 

More formally, we test the null hypothesis Hq : fj, — fio vs the alternate Hi : fj, ^ iiq ■ To do this, we 
proceed as follows: We compute the sample mean I and the sample standard deviation si. Then the standard 
deviation of I can be calculated as 

Sl 



s = 



Vn 



1 



n 

N 



(3) 



where a/1 — -^, the finite population correction term, is applied because we are drawing the subsample 
without replacement from a finite-sized population. Then, we compute the test statistic: 



I- P-o 



(4) 
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Figure 1: Theoretical (red line) and empirical distribution (blue bars) of the t-statistic under resampling n 
datapoints without replacement from a dataset composed of digits 7 and 9 from the MNIST dataset (total 
N = 12214 points). 

If n is large enough for the central limit theorem to hold, the test statistic t follows a standard Student- 
t distribution with n — 1 degrees of freedom, if the null hypothesis is true (see Figure |4] for empirical 
verification). Then, we compute the p- value S — 1 — 0„_i(|i|) where 0„_i(.) is the cdf of the standard 
Student-t distribution with n — 1 degrees of freedom. If 5 < e we can reject the null hypothesis, i.e. we have 
enough precision in the sample to make a decision. In this case, if Z > /io, we accept the proposal, otherwise 
we reject it. If we fail to reject the null hypothesis, we draw more data to reduce the uncertainty, s, in the 
sample mean I. We keep drawing more data until we have enough precision to make a decision. Note, that 
this procedure will terminate because when we have used all the available data, i.e. n = N the standard 
deviation s is 0, and the sample mean / = /i. So, in this case we will make the same decision as the original 
MH test would make. Pseudo code for our test is shown in Algorithm „ (here we have assumed that the total 



number of datapoints is divisi ble by the size of the mini-batch for simplicity) . This algorithm is related to 



the Pocock sequential design |Pocockl llQTTj which allows us to set a value of e given a particular value of 
Type I error. 

The advantage of our method is that often we can make confident decisions with n < N datapoints and 
save on computation, although we introduce a small bias in the stationary distribution. But, we can use 
the computational time we save to draw more samples and reduce the variance. The bias-variance trade-off 
can be controlled by adjusting the knob e. When e is high, we make decisions without sufficient evidence 
and introduce a high bias. As e — ^ 0, we make more accurate decisions but are forced to examine more data 
which results in high variance. 

Approximate MH test 
Input: 9t, 0' , e, m, Xn 
Output: accept 
1: Initialize estimated means T^— and P ^^ 

Initialize n •(— 

Draw u ^ Uniform[0,l] 

Compute Mo ^ ^ log [w^gljfjgf}] , 
done <— false 
while not done do 

Draw mini-batch Xm of size m without replacement from Xn 

Set X^ i — Xi\[ \ X^n 

Update I and P using Xm- 

n ^ n + m 

Update estimated standard deviation s -i— ^1 — -^ ( ~ _ 1 



2: 
3: 

4: 

5: 
6: 
7: 
8: 
9: 
10 
11 



Update student-t statistic t i ^^ 

Update p- value (5^—1 — 0„_i(|i|) 
if (5 < e then 

accept ^ true if Z > /io and false otherwise 

done <— true 
end if 
end while 



We will now prove an upper bound in the error of the stationary distribution of our approximate algorithm. 
Denote the transition kernel of the Metropolis-Hastings algorithm as T, and the total variation distance 

between two distributions, P, and Q as: dy{P,Q) = jll^" Qlli- 

Let P(accept|0, 9') be actual acceptance probability of our algorithm. Given an upper bound on the error 
of this probability: A^ax = supg g, |P(accept|0, 0') — Pa{0,9')\, when T satisfies a contraction condition, 
we can obtain an upper bound on the discrepancy between the posterior distribution and the stationary 
distribution of our approximate Markov chain as follows. 

Theorem 1. For a Metropolis-Hastings algorithm with the transition matrix T satisfying the following 
contraction condition 

d,{PT ~ S) < r]d,{P ^ S) (5) 

for any probability distribution P, the distance between the stationary distribution and that of the approximate 
sampler S^ is upper bounded by 



dy{So,Se) < 

The proof is provided in the appendix. 
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5 Application To Gibbs Sampling 

The same sequential test can be applied to the Gibbs sampling algorithm on discrete models with binary 
variables. Consider a probabilistic model on D binary variables P{Xi, . . . , Xd)- At every iteration of Gibbs 
sampling, it updates one variable Xi using the following procedure: 

1. Compute the conditional probability: 

where X-i denotes all the other variables. 

2. Draw u ^ Uniform[0, 1]. If u < P{Xi = l|a;-i) set Xi — 1, otherwise set Xi = 0. 
The condition in step 2 is equivalent to check 

log-u ^ hgPjX, ^ l,a:-,) 



log(l-u) \ogPiX,=0,x^^) 

When the joint distribution is expensive to compute but can be represented as a product over multiple terms, 
P{X) = n„=i fniX), we can apply the sequential test in the previo 
sampling algorithm. In this case the variable /j,0 and // is given by 



P{X) = n„=i fniX), we can apply the sequential test in the previous section in order to speed up the Gibbs 



1 V^ 1 fn{Xi ~ l,X^i) 

Similar to the Metropolis-Hastings algorithm, given an upper bound in the error of the approximate 
conditional probability 

Ainax = max|P(Xi is assigned l\x-i) - P{Xi = l|a::_i)| 
we can prove the following theorem: 



Theorem 2. For a Gibbs sampler with a Dobrushin coefficient 77 < 1 fBremaud . \ 19991 . chap. 7.6.2], the 



distance between the stationary distribution and that of the approximate Gibbs sampler S^ is upper bounded 
by 

1 — 7y 
The proof is also provided in the appendix. 

6 Experiments 

In this section we apply our approximate test to the following algorithms: 1) a random walk MH, 2) SGLD 
and 3) a Gibbs sampler. 

6.1 Random Walk proposal 

We first test our method using a random walk proposal q{9'\9t) = A/'(6't,CT|.^y). Although the random walk 
proposal is not efficient, it is very useful for illustrating our algorithm because the random walk proposal 
contains no information about the target distribution, unlike Langevin and Hamiltonian methods. So, the 
responsibility of converging to the correct distribution and sampling from it lies solely with the MH test. 
Also note that since q is symmetric i.e. q{9t\9') = q{9'\9t), it does not appear in the MH test at all. 



The target distribution in this experiment was the posterior for a logistic regression model trained on the 
MNIST dataset for classification of digits 7 vs 9. The dataset consisted of 12214 datapoints and we reduced 
the dimensionality from 784 to 50 using PCA. We chose a Gaussian prior on the parameters. 

In figure [U we show the percentage of data used(blue line) and the percentage of wrong decisions (red 
line) as a function of e. The amount of data required drops off very fast as e is increased even though the 
error does not increase much. Note that the setting e = corresponds to the exact MH algorithm. 
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Figure 2: Percentage of data (blue line) used and 
Percentage of wrong accept-reject decisions (red 
line) as a function of e. The black line shows e 
for comparison. Results are for a random walk 
proposal on a logistic regression model trained 
on the MNIST dataset 



Figure 3: Average Risk in predictive mean of 
2037 test points vs wall clock time for different 
values of e. T denotes the total number of sam- 
ples collected in 400 sees. Results are for a ran- 
dom walk proposal on a logistic regression model 
trained on the MNIST dataset 
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Figure 4: Marginals Ox vs Ot for e = and e = 0.01 at different values of wall clock time. The actual 
number of samples T is shown as well. The blue curve shows the true marginal obtained using a long run 
of Hybrid Monte Carlo. The red curve shows marginals obtained by running the random walk algorithm 3 
times. Results are for a random walk proposal on a logistic regression model trained on the MNIST dataset 



In figure m we show marginals of 6i vs 62 for e = and e — 0.01 at different values of wall clock time. 
The red curves are marginals for 3 different runs of the random walk algorithm whereas the blue curve shows 
the true marginal obtained from a long run of Hybrid Monte Carlo. At T = 50 sec, the exact MH algorithm 



e = has still not completed burn-in. Our algorithm with e — 0.01 is able to accelerate burn-in because 
it can collect more samples in a given amount of time. Theoretically, as more computational time becomes 
available the exact MH algorithm will catch up with our algorithm, and eventually dominate it, because the 
exact MH algorithm is unbiased unlike ours. But note that the bias is hardly noticeable in this example and 
the error is still dominated by the variance even after collecting around lOOK samples. 

In figure [31 we show how the logarithm of the risk in the predictive mean of a test point decreases as a 
function of wall clock time. The risk is computed as the mean squared error using 8 different simulations of 
the Markov chain, and we average this risk in predictive mean over 2037 datapoints in the test-set. Since we 
are plotting log R = log{B^ + V) = \og{B^ + ^jt-), we expect to see a straight line with a negative slope when 
the variance dominates the bias. Eventually as the bias starts to dominate we would expect slope of these 
lines to go to zero. The figure shows that even after collecting a lot of samples, the risk is still dominated 
by the variance. The minimum risk after a finite amount of computational time is obtained with e > 0. 

6.2 SGLD proposal 

We also combined our method with the SGLD kernel. The SGLD proposal is: 



2^^ 



J2 \ogpix\9) + pi9) 



(10) 



.xex„ 



where i] ^ Af(0,a). and Xn is a random mini-batch of size n. Since the proposal depends on a randomly 
drawn mini-batch, we can think of SGLD as being a mixture kernel, where each mini-batch corresponds to 
a different kernel. To ensure that a mixture kernel satisfies detailed balance, we just have to ensure that 
each each kernel satisfies detailed balance. Thus, in the MH formula we use q{6t\9',Xn) and q{d'\6t,Xn) 
conditioned on the same mini-batch. 

In this experiment, we trained a logistic regression classifier on the KDD99 dataset to detect malicious 
network intrusions. The dataset contains about 4 million datapoints, but we train on a 10 % subset. There 
are 41 parameters in the model. 
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Figure 5: LI error in variance vs Wall Clock Time 
with SGLD proposal. Step size = 5e-4 



Figure 6: LI error in variance vs Wall Clock Time 
with SGLD proposal. Step size = le-4 



We compare uncorrected SGLD, MH corrected SGLD and SGLD corrected using our sequential test in 
Figures [5] and [6] . We plot the LI error in variance as a function of wall clock time at two different step sizes 
for the SGLD kernel. The results are averaged over 8 Markov chain runs. Note that when the step size is 
large SGLD has a huge bias. While the MH test can correct this, the MH test introduces a huge variance. 
Our method on the other hand can reduce the bias while at the same time maintaining a low variance. On 
the other hand when the step size is small, the bias in SGLD is small. Since SGLD is very fast, it can reduce 
the error quicker than other methods. But, if we run this longer, the other methods can catch up and achieve 



a lower error than SGLD. Note however that for both figures, the MH test is prohibitively slow, although 
given an infinite amount of time it can achieve error while other methods cannot. 

These results suggest that the ideal case would be a combination of these methods. We can start 
with uncorrected SGLD to achieve fast convergence to a high probability region. Then we can add our 
approximate MH tests to ensure that detailed balance is approximately satisfied, while still not being too 
much of a computational burden. As more samples are collected, the variance reduces. So, we can keep 
reducing e to balance the contributions of bias and variance to the overall error. Eventually in the limit of 
infinite computational time, we approach the true MH algorithm. 



6.3 Gibbs Sampling on Markov Random Fields 

We illustrate the performance of our approximate Gibbs sampling algorithm on a synthetic Markov Random 
Field. The model under consideration has D = 100 binary variables and they are densely connected by 
potential functions of three variables ^ij^k{Xi,Xj,Xk),yi y^ j y^ k. There are totally D{D — 1){D — 2)/6 
potential functions (we assume potential functions with permuted indices in the argument are the same 
potential function), and every function has 2"^ = 8 values. The entries in the potential function tables are 
drawn randomly from a log-normal distribution, log ipij^k{Xi,Xj,Xk) ^ A/'(0, 0.02). To draw a Gibbs sample 
for one variable Xi we have to compute (D — 1){D — 2)/2 = 4851 pairs of potential functions as 



p{x, = i|x_,) n 



i^j^k 



1pi,].kiXi = l,Xj,Xk) 



P{Xi = 0\x-i) lli^-i^k i^i,j,kiXi = 0,Xj,Xk 



(11) 



The approximate methods use a mini-batches of 500 pairs of potential functions at a time. We compare the 
exact Gibbs sampling algorithm with approximate versions with e G {0.01, 0.05, 0.1, 0.15, 0.2, 0.25}. 

To measure the performance in approximating P{X) with samples Xt, the ideal metric would be a distance 
between the empirical joint distribution and P. Since it is impossible to store all the 2^°° probabilities, we 
instead repeatedly draw M = 1600 subsets of 5 variables, {sm}m=i' s™ C {1, . . . , D}, \sm\ = 5, and compute 
the average Li distance of the joint distribution on these subsets between the empirical distribution and P: 



Error=^5]||P(X3„)-P(X,„) 
The accurate P is estimated by running Gibbs chains for a long time. 



(12) 



0.8 



5 0.6 
o 

qI 
•^0.4 

Q. 

E 

LU 





jf- 


^ 


^^^ 


True probability 


y^^ 


e = 0.01 


/^;:f0^ 


E = 0.05 


/^^^^^^^^"^ 


E = 0.1 


/0^ 


E = 0.15 

E = 0.2 


r 


E = 0.25 



1.4 

1.2 

1 

o 0.8 

L 0.6 

0.4 

0.2 



0.2 



0.2 0.4 0.6 0.8 1 

True Conditional Probability 

Figure 7: Empirical conditional probability vs exact 
conditional probability for different values of e. The 
dotted black line shows the result for exact Gibbs sam- 
pling. 
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Figure 8: Average L\ error in the joint distribu- 
tion over cliques of 5 variables vs running time 
for different values of e. The black line shows the 
error of Gibbs sampler with an exact acceptance 
probability. T in the legend indicates the number 
of samples obtained after 1000 seconds. 



We show the empirical conditional probability obtained by our approximate algorithms (percentage of 
Xi being assigned 1) for different e in Figure [T] It tends to under estimate at large probabilities and over 
estimate on the other end. When e = 0.01, the observed maximum error is within 0.01. 

Figure |5] shows the error for different e as a function of the running time. For small e, we use fewer 
mini-batches per iteration and thus generate more samples in the same amount of time than the exact Gibbs 
sampler. So the error decays faster in the beginning. As more samples are collected the variance is reduced. 
We see that these plots converge towards their bias floor while the exact Gibbs sampler out-performs all the 
approximate methods at around 1000 seconds. 

7 Conclusions and Future Work 

In this work we have taken a first step towards cutting the computational budget of the Metropolis-Hastings 
MCMC algorithm, which takes 0{N) likelihood evaluations to make the binary decision of accepting or 
rejecting a proposed sample. In our approach we compute the probability that a new sample will be accepted 
based on a subset of the data. We increase the cardinality of the subset until a prescribed confidence level 
is reached. In the process we create a bias, which is more then compensated for by a reduction in variance 
due to the fact that we can draw more samples per unit time. Current MCMC procedures do not take these 
trade-offs into account. In this work we have focused on a fixed decision threshold for accepting or rejecting 
a sample, but in theory a better algorithm can be obtained by adapting the decision threshold over time. 
An adaptive algorithm can tune bias and variance contributions in such a way that at every moment our 
risk (the sum of squared bias and variance) is as low as possible. We leave these extensions for future work. 

A Proofs of Theorem 1 and 2 

A.l Notations 

Denote P{X) as a probability distribution on the state space X. Let T be the transition kernel of a Markov 
chain, and % be the kernel of the approximate sampling algorithm. When we simulate a Markov chain from 

an initial distribution P, denote the distribution after t steps by P*-*' = PT^*^ 

A. 2 A Lemma for Theorem 1 and 2 

We first prove a lemma that would be used for the proof in both theorem 1 and 2. 

Lemma 3. Given two transition kernels, T and T, with respective stationary distributions, S and S , if T 
satisfies the contraction condition with a constant rj G (0, 1) as follows: 

d,iPT-S)<T]d^{P^S),yP (13) 

and the one step error between T and T is upper bounded uniformly with a constant S > as 

d,iPT - Pf) < 6,yP (14) 

then the distance between S and T is bounded by 

d.(5-5)<-i- (15) 

l-?7 

Proof. Consider a Markov chain with the transition kernel T, started from an arbitrary distribution P: 
p(o)^ p(i)^ .... At any time step, t > 0, we apply the transition kernel T on P^*\ According to the one step 
error bound, the distance between p(*+i) and the distribution obtained with T is upper bounded by 

d„(p(*+i)-p(*)r)<(5 (16) 

Following the contraction condition of T, the distance of P'^^^T from its stationary distribution 5* is less than 
P(*) as 

d„(p(*)r- 5) < ryd.(P(*' - S) (17) 

10 



Now let us use the triangular inequality to combine Equation 1161 and 1171 to obtain an upper bounded for the 
distance between p(*+i) and S: 

d,(p(*+i) -S)< rf,„(p(*+i) - P^'^T) + d,{P'^''>T - S) 

<<5 + 77d„(p(*)-5) (18) 

For any positive constant r < 1 — 77, when dy{P^*'^ — S) > ^, <^ 5 < rdy{P^*^ — S), plug it into Equation [TBI 
and we get a contraction condition for T towards 5: 

d^(p(*+i) -S)<{r + rj)dy{P^*') - S) (19) 

So if the initial distribution P is outside the ball B{S, -) == {p : dy{p, S) < -} in the space of distributions, 
the Markov chain will move monotonically into the ball at a finite number of steps for any r < 1 — 77. 
Denoting the first time it enters the ball as T(r). If P is already in ^(5*, -), we simply let T{r) — 0. We then 
show by induction that P*^*-* will stay the ball for all t > T{r). 

1. At t = T{r), P(*) e 8(8, ^) holds by the definition of r. 



2. Assume P^*' G B{S, -) holds for some t > T{r). Then following Equation [181 we have 



d„(p(*+i)-5') <<5 + ry- 



r + 11 6 

-S < - 

r r 



So P'-*-' e 6(5', -) holds for any t > T(r). Since P^*) will converge to S, we know that 

dy{S,S) < -,Vr <l-ry 
r 

Taking the limit r — >■ 1 — 77, we prove the lemma as 

dyiS,S)< 



1 — r] 



(20) 

(21) 

(22) 
D 



A. 3 Proof of Theorem 1 

We first derive an upper bound for the one step error for the Metropolis-Hastings algorithm, and then use 
Lemma [3| to prove Theorem 1. 

The transition kernel of the Metropolis-Hastings algorithm can be written as 

r(x, y) = Pa{x, y)q{y\x) + (1 - Pa{x, y))5D{y - x) (23) 

where 5d is the Dirac delta function. For the approximate algorithm proposed in this paper, we have an 

approximate acceptance function Pa(a;, y) where the error, APa = P—P, is upper bounded by | AP^I < A,„ax- 

Now let us look at the distance of the distributions generated by one step of T and the approximate 

kernel T. For any P, 

WPf-PTh 

dPix)APa{x, y) {q{y\x) - Soiy - x)) 

1 

< A„ 



— 2A 



dPix){q{y{\x) + Soiy - x)) 

Q{y) + P{y)\\i 



def 



where Q{y) — J dP{x)q{y\x). So we get an upper bound for the total variation distance as 

dy{Pf,PT) = \\\Pf-PT\\i < A^ax 
Apply Lemma [3] with 6 = A^ax and we prove Theorem 1. 



(24) 



(25) 



11 



A. 4 Proof of Theorem 2 

The proof of Theorem 2 is similar to that of Theorem 1. We first obtain an upper bomid for the one step 
error and then plug it into Lemma [31 

The transition kernel of the Gibbs sampler for variable Xi can be represented by a matrix 71 of size 
2^ X 2^: 

n{x,y) - I p^Y^ ^ ^^|^_^^ otherwise ^^> 

where 1 < i < N,x,y Q {0, 1}^. The approximate transition kernel T can be represented similarly as 

'^^ ' ^^ " I P{Y. = y^\y-^) otherwise ^^'> 

where P is the approximate conditional distribution. Define the approximation error /S.T{x, y) = T{x, y) — 
T{x,y). We know that AT{x,y) = if y_, ^ X-i and it is upper bounded by Amax from the condition of 
Theorem 2. 

For any distribution P, the one step error is computed as 

\\Pf-PT\\i 



E 

V 

E 



^P(x)Ar(a;,j/) 



2:,e{0,l} 



So we get an upper bound for the total variation distance 



(28) 



d,{Pf,PT) - \\\Pf-PT\\i < A„,ax (29) 



For a Gibbs sampling algorithm, we have the contraction condition [Bremaudl Il999l chap. 7.6.2] as 



d,iPT - S) < r]d,{P ~ S) (30) 

where rj is the Dobrushin ergodic coefficient. Plug S — Amax and rj into Lemma [3] and we obtain the 
conclusion in Theorem 2. 

References 

S. Ahn, A. Korattikara, and M. Welling. Bayesian posterior sampling via stochastic gradient fisher scoring. 
In International Conference on Machine Learning, 2012. 

C. Andrieu and G. O. Roberts. The pseudo-marginal approach for efficient monte carlo computations. The 
Annals of Statistics, 37(2):697-725, 2009. 

P. Brcmaud. Markov chains: Gibbs fields, Monte Carlo simulation, and queues, volume 31. Springer, 1999. 

D. Gamerman and H. F. Lopes. Markov chain Monte Carlo: stochastic simulation for Bayesian inference, 
volume 68. Chapman & Hall/CRC, 2006. 

Z. Govindarajulu. Sequential statistics. World Scientific Publishing Company Incorporated, 2004. 



12 



L. Lin, K. Liu, and J. Sloan. A noisy monte carlo algorithm. Physical Review D, 61(7):074505, 2000. 

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of state calcula- 
tions by fast computing machines. The journal of chemical physics, 21:1087, 1953. 

S. J. Pocock. Group sequential methods in the design and analysis of clinical trials. Biometrika, 64(2): 
191-199, 1977. 

S. Singh, M. Wick, and A. McCallum. Monte carlo mcmc: efficient inference by approximate sampling. 
In Proceedings of the 2012 Joint Conference on Empirical Methods in Natural Language Processing and 
Computational Natural Language Learning, pages 1104-1113. Association for Computational Linguistics, 
2012. 

M. Welling and Y. Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 
28th International Conference on Machine Learning (ICML), pages 681-688, 2011. 



13 



